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Time reversal of acoustic waves can be achieved efficiently by the persistent control of excitations 
in a finite region of the system. The procedure, called Time Reversal Mirror, is stable against the 
inhomogeneities of the medium and it has numerous applications in medical physics, oceanography 
and communications. As a first step in the study of this robustness, we apply the Perfect Inverse 
Filter procedure that accounts for the memory effects of the system. In the numerical evaluation of 
such procedures we developed the Pair Partitioning method for a system of coupled oscillators. The 
s ! \ algorithm, inspired in the Trotter strategy for quantum dynamics, obtains the dynamic for a chain 

of coupled harmonic oscillators by the separation of the system in pairs and applying a stroboscopic 
' • sequence that alternates the evolution of each pair. We analyze here the formal basis of the method 

fj^ | and discuss his extension for including energy dissipation inside the medium. 



I. INTRODUCTION 

13 . 

In the recent years, the group of M. Fink developed an experimental technique called Time Reversal Mirror (TRM) 
[l[ that allows time reversal of ultrasonic waves. An ultrasonic pulse is emitted from the inside of the control region 
(called cavity) and the excitation is detected as it escapes through a set of transducers placed at the boundaries. 
These transducers can act alternatively like microphones or loudspeakers and the registered signal is played back in 
the time reversed sequence. Thus, the signal focalizes in space and time in the source point forming a Loschmidt Echo 
0. It is remarkable that the quality of localization gets better by increasing the inhomogeneities inside the cavity. 
This property allows for applications in many fields [H, 0] • In order to have a first formal description for the exact 
reversion, we introduced a time reversal procedure denoted Perfect Inverse Filter (PIF) in the quantum domain []|. 
The PIF is based in the injection of a wave function that precisely compensates the feedback effects by means of the 
renormalization of the registered signal in the frequency domain. This also accounts for the correlations between the 
q ' transducers. Recently, we proved that these concepts apply for classical waves @. We applied it to the numerical 
evaluation of the reversal of excitations in a linear chain of classical coupled harmonic oscillators with satisfactory 
results. A key issue in assessing the stability of the reversal procedure is have a numerical integrator that is stable and 
perfectly reversible. Therefore, we developed a numerical algorithm, the Pair Partitioning method (PP), that allowed 
the precise test of the reversal procedure. In the next section we develop the main idea of the method and then 
\ we use it to obtain the numerical results that we compare with the analytical solution for the homogeneous system. 
y—i . Additionally, we introduce a method to approximate the solution of an infinite system using a finite one. For this we 
introduce a non-homogeneous fictitious friction term that can simulate the diffusion of the excitation occurring in an 
^vq , unbounded system. These strategies are tested through a numerical simulation of a time reversal experiment. 

o : 

00 . 

O ■ II. WAVE DYNAMICS IN THE PAIR PARTITIONING METHOD 

> : 

The system to be used is shown in the figure Q] a one-dimensional chain of N = 2s coupled oscillators with masses 
rS [ rrii and natural frequencies u>i that can be represented as a set of coupled pendulums. 

If pi denotes the impulse and Ui the displacement amplitude from the equilibrium position for the ith oscillator, 
the Hamiltonian writes 



> 
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where Ki is the elastic coefficient that accounts for the coupling between the oscillators i and i + Notice that we 
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FIG. 1: Scheme of the N coupled oscillators system to be solved by the PP numerical method. 



could rewrite the Hamiltonian in terms of each coupling separating it in non-interacting terms including even pairs 
and odd pairs each: 



T~L — ~Hodd + Woven — ^1,2 + ^3,4 + • • • + Hn-1,N 
+ Ti-2,3 + 71-4,5 + . • • + Hn-2,N-1, 
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and 



m 1 =m u lo 1 =lo 1 , m N = m N , u N =lo n , 

(4) 

m n = 2m„, w„ = w n /2, n = 2, . . . , Jv — 1. 

A good approximation to the overall dynamics, inspired in the Trotter method used in quantum mechanics 0], 
can be obtained solving analytically the equations of motion for each independent Hamiltonian in a time step r. 
Therefore, the pair 7i nin+ i has 



- 2 , \ K n 

u n + ^— u n + — U n+ i, 

m n I m„ . . 

7 . (5) 
, K n \ K n 

in+l = - U n+1 + — U n+1 + U n . 

\ m n +ij m n+ i 

At each small time step r, the evolution for the even couplings is obtained and the resulting positions and velocities 
are used as initial conditions for the set of Hamiltonians accounting for odd couplings and so on. Since the equations 
of motion are solved separately, we could consider only the two coupled oscillators system, e.g. 



u 2 = -u\u 2 + Ul\ x U\. 
For this system, it is easy to obtain the corresponding normal modes 



(6) 
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u ± = T-2 — -r«i ± 7-5 — -r u 2, (7) 



with characteristic frequencies 



From ||7J|, U±(t) are obtained and these values are used for the evolution after the temporal step t — > t + r, i.e. 

0±(* + r) = J7 ± (t)cos(w ± r) + L>±(t)sin(w ± r)/w ± , 

(9 J 

f/±(i + r) = L/±(f) cos(w±t) — U±(t) 8w{w±t)w±. 
Once we have U±(t + r) and U±(t + r) we can go back to the natural basis by means of the inverse of j7]) 

22 22 

UJ\ — 0J± tof — OJ~L 

ui = * U+ H n £/_, . . 

Wll ^1 (10) 

1X2 = £/+ + C/_. 

Then, one obtains the displacements and momenta for all oscillators at time i + r. The above steps are summarized 
in the Pair Partitioning algorithm: 

1. Determine all the masses and natural frequencies of the partitioned system ifi n ,uj n . 

2. For even couplings 7^2,3,7^4,5, ■ ■ •, rewrite the initial conditions {ttj(0), ttj(Q)} for the normal modes according 
to 0. 

3. Calculate the normal modes evolution for even couplings, according to (0 and obtain {u° von (r), u° von (r)} from 



4. Calculate the normal modes for odd couplings Hi, 2, 7^3.4, . . ., using the recent positions and velocities. 

5. Calculate the normal modes evolution for odd couplings and give the positions and velocities {ui(r), Ui(r)}. 

6. Go back to the step 2 with {u,-(t),u,(t)}. 

Therefore, applying n times the PP algorithm we obtain the positions and velocities for all oscillators at time 
t = nr. 

As an example we consider the homogeneous system where the N oscillators have identical masses and the only natural 
frequencies correspond to the surface u>i = lon = uj x = K/m. The displacement amplitude of the ith oscillator due 
to an initial displacement in the jth oscillator can be expressed analytically as 

uftS) = — - — y^sin ( i — — ) sin (j — - — ) cos(w fc *), (11) 

N+l^ \N + lJ \ J N+lJ V h y ' 

with 

u>k = w x sin ( — -1 , (12) 

\2(N + 1))' V ' 

the characteristic frequency for the fcth normal mode. In the figure [2] the analytical and numerical results are 
compared for the surface oscillator displacement in a case when all the oscillators were initially in their equilibrium 
positions except ui(0) = uj. We use N = 200 in two cases where the temporal steps are r = lO -2 ^" 1 and 

r = 10-Vr 1 . 
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FIG. 2: Time evolution for the surface oscillator. Top, comparison between the analytical and numerical results shows all 
curves superposed. Bottom, error in the displacement for both temporal steps. 

We have taken N such that no mesoscopics echoes @ appear in the interval of time shown. We also notice that 
the error 



e (t) = 1 1 w 13 

drops as the temporal step r diminishes. We observe a quadratic dependence maxe(i) = ar 2 in complete analogy 
with the Trotter method. In the particular case of the homogeneous system we have a ~ 0.0445. 



A. Unbounded systems as damped oscillations 



The solution of wave dynamics in infinite media remains as a delicate problem. In such case, the initially localized 
excitation spreads through the systems in a way that resembles actual dissipation. In contrast, finite systems present 
periodic revivals, the mesoscopic echoes, that show that energy remains in the system. In order to get rid the 
mesoscopic echoes and obtaining a form of "dissipation" using a finite number of oscillators, we add a fictitious 
"friction" term. The friction coefficients rji > can be included between the 2th and 3th steps of the PP algorithm 
supposing that the displacement amplitude decays exponentially 



u,i{t) -> Ui(t) exp(-^r), (14) 

as occurs in a damped oscillator in the limit rj/uj <C 1. For the homogeneous system with N = 100 oscillators were 
the cavity ending at site xr we choose a progressive increase in the damping as 

th = 0.1 * ~ Xa , i = x R ,...,N. (15) 
N - x R 

We compare the result of this approximation with the undamped case for the displacement amplitude in xr = 10. 
The figure [3] shows how the dynamics in the damped system has no mesoscopic echoes whereas in the undamped 
system we observe the echo at t M ~ 2Ntu~ 1 . 

As we will see for TRM and PIF procedures, this last result is very usefull since it allows to obtain the dynamics 
of an open system with a small number of oscillators (e.g. N ~ 100). 
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FIG. 3: Comparison between the undamped evolution (black) and the damped evolution (blue) for the displacement amplitude 

Wil(t). 

III. NUMERICAL TEST FOR THE PERFECT INVERSE FILTER PROCEDURE 

As we mentioned above, the Time Reversal Mirror procedure consists in the injection, at the boundaries of the 
cavity, of a signal proportional to that recorded during the forward propagation. In contrast, the Perfect Inverse 
Filter corrects this recorded signal to accounts for the contributions of multiple reflections and normal dispersion 
of the previously injected signal in a manner that their instantaneous total sum coincides precisely with the time 
reversed signal at the boundaries. The continuity of the wave equation ensures that perfect time reversal occurs at 
every point inside the cavity. The procedure for such correction is described somewhere else [H, H| . Here, it is enough 
to notice that the imposition of an appropriate wave amplitude at the boundaries should give the perfect reversal of an 
excitation originally localized inside the cavity. An example of such situation would be a " surface pendulum" coupled 
to a semi-infinite linear chain of harmonically coupled masses |9| . In such a case, we know that the energy decays in 
a approximately exponential way, where the decay rate can be assimilated to a "friction coefficient " . However, for 
very short times, the local energy decays with a quadratic law while for very long times the exponential decay gives 
rise to a power law characteristic of the slow diffusion of the energy [l(| • In the figure 0] we show how this overall 
decay is reversed by controling the amplitude at site xr = 10. 

As long as we have been able to wait until a neglegible amount of energy is left in the cavity, the control of the 
boundaries is enough to reverse the whole dynamics inside the cavity. As a comparison, the theoretical reversal of the 
decay of a surface excitation in a semi-infinite chain is shown. We see that the region they differ is when the injected 
signal is still negligible as compared to the energy still remaining in the cavity (as a concequence of the very slow 
power law decay). 

IV. DISCUSSION 

We have presented a numerical strategy for the solution of the wave equation that is completely time reversible. 
This involves the iterative application of the exact evolution of pairs of coupled effective oscillators where the energy 
is conserved, hence deserving the name of Pair Partitioning method. This is complemented with an original strategy 
for dealing with wave propagation through infinite media. While various tests of these procedures remain to be done, 
we have shown, through the solution of simple but highly non-trivial examples, that the method is numerically stable 
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FIG. 4: Local energy of the surface oscillator recovering. Here, we choose for the detection time tn — lOOOr. evolution (blue) 
for the displacement amplitude un(t). 



and can be used to revert wave dynamics up to a desired precision. 
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